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We develop a three-dimensional analysis of the phase sep- 
aration of two-species Bose-Einstein condensates in the pres- 
ence of vorticity within the Thomas-Fermi approximation. 
We find different segregation features according to whether 
the more repulsive component is in a vortex or in a vortex-free 
state. An application of this study is aimed at describing sys- 
tems formed by two almost immiscible species of rubidium-87 
that are commonly used in Bose-Einstein condensation exper- 
iments. In particular, in this work we calculate the density 
profiles of condensates for the same conditions as the states 
prepared in the experiments performed at JILA [ Matthews, 
et ai, Phys. Rev. Lett. 83, 2498 (1999)]. 
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Vortices in Bose-Einstein condensates (BECs) have 
thoroughly been studied from a theoretical point of view 
since the first experiments on condensates were carried 
out, due to their direct connection to superfluidity. For a 
review of these issues see for example Q , and for the most 
recent works references in Q]. But until the experiment 
of Matthews et al. [|j was performed, no evidence about 
their existence had ever been reported. This experiment 
was based on a method proposed by Williams and Hol- 
land 1^ to create vortices in two-component BECs. They 
trapped atoms of ^^Rb in two different hyperfine states 
\F — l,mf = —1) and \F = 2,m/ = 1), henceforth 
denoted as |1) and |2), respectively, and created vortices 
either in the |1) or |2) component, while the other species 
remained without rotation. Because of the relation be- 
tween intra- and inter-particle scattering lengths, these 
states have the property of being almost immiscible. We 
shall see that this fact has important consequences on 
the structure of the vortices. 

In connection with this experiment interesting theo- 
retical work has recently been developed in order 
to describe different properties of many types of two- 
component condensate systems. Garci'a-Ripoll and Perez 
Garcia, by means of a robust formalism, analyzed the sta- 
bility and dynamics of a variety of two-species conden- 
sates, taking into account the effect of the interchange of 
species in the system, and also the variation of the num- 
ber of particles. In particular, in Ref. |^ they set the 
relative number of particles equal in both species and 
varied the total number of particles, and in Ref. |^ they 
pursued the research by also varying the relative compo- 
sition using a 2D model. However, in these works [||J|] 



little is said about the structure of the stationary states 
and no comments are made about segregation of species. 
In a recent work, Chui et al. computed the energy 
of different configurations within the Thomas- Fermi Ap- 
proximation (TFA). They varied the number of particles 
and interchanged the species in order to derive conclu- 
sions about the relative stability. But, although they 
used a 3D model and considered phase separation in the 
mixture, a hypothesis was made regarding the form of the 
interface between species that we find is not consistent 
with our results and apparently also with new available 
experimental data We shall comment on this later. 

In 1996, in a pioneering work, Ho and Shenoy Q con- 
structed an algorithm to determine the density profiles 
of binary mixtures of alkali atoms within the TFA. We 
develop here a different procedure which allows one to 
obtain conclusions about the general features of segrega- 
tion through simple calculations. 

The aim of the present article is to analyze and per- 
form a classification of the structure of condensates near 
phase separation, when two almost immiscible species are 
involved. Our starting point for calculating the station- 
ary states are the Gross-Pitaevskii equations in the TFA 

0: 

= (14+iViGi,i|*i|2+iV2Gi,2|*2n*l, 
/i2*2 = {V2 + iViGi,2|*lP + iV2G2,2|*2n *2, (1) 

where Ni denotes the number of atoms of species i, and 
Gk,i — Uk,iU, with Uk^i the relative interaction strengths 
between species k and I. We set the most repulsive com- 
ponent |1) in the '^i state, and fix ui i = 1 (> U2,2), 
so U = -^-5(7^1 being the scattering length of the |1)- 
species and M the atom mass. 

For simplicity, we consider the system in a spherically 
symmetric trap and make a change of variables accord- 
ing to \J~^-wof f, with Wo the trap angular frequency. 
Thus the effective potentials Vi written in cylindrical 
variables read 

V, = r^ + z^ + ^. (2) 

The first two terms correspond to the trapping poten- 
tial, while the third includes the centrifugal term 

Ki = ™° , with nii the number of quanta of circu- 
lation associated with the wave function for species i. 
Let us first analyze what happens when D = U1.1M2.2 — 
2 = 0- Physically the system starts undergoing phase 
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separation. Mathematically, for both > 0, Eqs. 
become linearly dependent. So to have a compatible set 
of equations the following condition must hold: 
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If we define P2 — '"i,2/'"i,ii the above equation can be 
rewritten as 



H2 - P2H1 
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This equation defines different families of surfaces that 
we shall call Sq, where i labels the family. As was men- 
tioned above, this is the only region of space where both 
^'i may be non-zero, so these surfaces necessarily deter- 
mine the interfaces between species. Such families are 
characterized by the values of 
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(note that D = implies P2 < 1)- 

Depending on the values of the above parameters, 
three topologically different kinds of surfaces exist, 
namely, spheres, cylinder-like surfaces, and toroids. In 
Fig. 1 we show a partition of the parameter space {A, B) 
indicating the domains corresponding to each family. In 
order to classify these surfaces let us study the z(r) curves 
defined by Eq. (^). If z{r) intersects the z axis, the asso- 
ciated surfaces 5g are spheres, and this occurs for i? = 
and A > Q. If z{r) has only one root the corresponding 
surfaces 5q are cylinder-like, and this takes place when 
B > Q. Finally, if z{r) has two roots the related surfaces 
Sq are toroids, and it is easy to verify that the conditions 
to be fulfilled are -t- 4B > 0, A > 0, and B <Q. The 
remaining region of the parameter space has no associ- 
ated surfaces. 

To perform a displacement of a point within the pa- 
rameter space two mechanisms are possible. On the one 
hand, the A value may be changed by varying the number 
of particles in each species, which will cause a change in 
the chemical potentials. A variation of this type does not 
change the family, but the shape of the surface is mod- 
ified. For example, for _B > one can get a straighter 
"cylinder" by moving from positive to negative A values, 
and, in terms of the number of particles, this corresponds 
to a decrease in the N2/N1 ratio. On the other hand, B 
can only take discrete values and it is modified by chang- 
ing the vorticity. 

It is worthwhile mentioning that the only possible real- 
ization for _B = 0, (and hence for 5q), within reasonable 
values of m^, is a vorticity free condensate mi = m2 — 0. 
In such a case, regarding the location of each species, it 
is easy to conclude that the energy minimum is obtained 
when the less repulsive component lies inside the sphere 
of radius VA, while the other species forms a spherical 



shell around it. In the presence of vorticity (_B ^ 0), the 
location of each species can be predicted from energetic 
considerations. The state with greater vorticity should 
be located in the region of space that excludes r = 0, 
because of the centrifugal term. 

From this simple analysis we may state that close to 
phase separation the presence of vorticity dramatically 
changes the segregation structure. 

Now, let us consider the general case £> 7^ 0. The solu- 
tion of Eqs. (0 ) can be easily obtained and has the fol- 
lowing different expressions depending on whether there 
exists any overlap between the wave functions of both 
species, a) In the region where only one wave function is 
non-vanishing ( ^ and l^'fcp = 0, for i ^ k) Eqs. 

(|l|) are decoupled and the solution is 
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b) In the region where both wave functions are non- 
vanishing, l^'ip > 0, one obtains 
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where /3i = ^ and A, - 

We call Si { S2) the surface defined by equating the 
expression inside the square bracket in |\I'ip (|^'2p) to 
zero. It is interesting to notice that the expressions men- 
tioned above are of the same type as the one that defines 
So (cf. Eq. (|4|)), and so a similar classification can be 
applied to these surfaces. The surfaces Si and ^2 define 
the boundary of the coexisting region of the |1) and |2) 
species. Depending on the value of the interparticle in- 
teraction strength, the regions where the overlapping of 
the waves function occurs exhibit different features. For 
D > and ui,2 small enough, the coexisting region can 
be very large, as may be observed in the figures shown in 
Rcf. 0. When increasing ui,2 this region becomes nar- 
rower. At = 0, it is easy to verify that 51=52 = 
Sq , and hence the coexisting region is reduced to a single 
surface. When ui 2 increases making D < 0, this surface 
splits into two distinct surfaces, forming thus a new nar- 
row 3D overlapping region, due to mutual repulsion of 
the species. 

To have the problem completely solved one has to cal- 
culate the chemical potentials. In order to obtain /ii and 
fi2 the normalization condition for the wave functions 
given by Eqs. (^) and (^ must be used. This leads to 
a system of two coupled integral equations that must be 
solved numerically. We have accomplished this in order 
to describe the states obtained in the JILA experiment 
Q, choosing parameters equal to those of the experi- 
ment. Namely, the trapping potential has an angular 
frequency wq = 27r ftrap, being ftrap = 7.8 Hz. The inter- 
action strengths relative to the |1) component of the ®^Rb 
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species are iti^ = l,ui,2 = -97 and U2.2 = -94. And the 
number of particles of each species is Ni = N2 = 4x10^. 
Note that for these values of the interaction strengths, 
D = —9x10"'' is close to the phase separation value, so 
the region between interfaces is very narrow, and hence 
Si ^ S2 ^ Sq. Thus hereafter, in order to identify the 
type of segregation of species encountered, we will make 
use of the classification we performed for Sq- 

For completeness, we also include calculations for the 
two species without vorticity, although the experimen- 
tal realization consisted only of one vortex and one non- 
vortex state in each component. 

In Fig. 2 we show the density contours of each wave 
function in the (r,z) plane for the following three con- 
figurations, (i) A vortex-free system (rni,m2) = (0,0) 
yields chemical potentials in units of fiwo as /ii — 24.50 
and /i2 = 24.06. It may be seen that in this case the 
segregation evolves as ~ Sq. The more repulsive compo- 
nent is on a spherical shell around the other, as expected, 
(ii) For (mi,TO2) = (1,0) we obtained /ii — 24.54 and 
fi2 — 24.04. The vortex segregates in a ring surrounding 
the non vortex state. It is easy to verify that the species 
become separated by an interface of the type Sq. The 
core size is about 9/xm and the whole vortex has a radius 
of 26/xm. This reproduces very well the experimental 
data, (iii) When the vortex is in the less repulsive com- 
ponent (mi,m2) = (0,1), we obtained /ii — 24.50 and 
/j,2 = 24.15. In this case the vortex is segregated inside 
the non- vortex state. The corresponding surface that de- 
scribes such segregation is the toroid Sq. The inner and 
outer radii of the torus are 5/xm and 18/im, respectively. 
However, this does not correspond to what was obtained 
at the first stage of the experiment, which looks quite 
similar to the pattern described in (ii). Mainly, in the 
experiment, the radius of the vortex was approximately 
26/zm, while according to our calculations it should be 
18/im. Taking into account that after vortex creation, 
the vortex immediately shrank, a possible and likely ex- 
planation for the discrepancy is that the system was not 
created near the equilibrium configuration but moved to- 
wards it. The case with a vorticity of (mi, 7712) = (1,0) 
is created much closer to equilibrium |lO[ . It is impor- 
tant to notice that in every case, even when the vortex 
configuration is out of equilibrium, the whole condensate 
is almost spherical with the same radius value of 26/im. 

In Fig. 3 we display the density profiles as a function 
of the radial coordinate at z = 0, for the same three 
configurations described in Fig. 2. It may be seen in 
all of the cases that the shape of the total density is 
like a parabola, a characteristic that also holds for any 
other direction considered. Regarding the densities of 
each one of the components separately, they abruptly 
go to zero in the narrow region between species. It is 
worthwhile to mention that this behavior would not be so 
sharp for real systems and some interpenetrating region 
would exist 1 1 1 1 . 



With respect to the hypothesis made by Chui et al. 
, they assumed for condensates including vorticity in a 
symmetric trap that one component is confined inside a 
ball and the other one in a spherical shell around it; and 
that the species with vorticity have a small healing length 
sized vortex core. In the context of our article this means 
that they worked as if segregation of species, including 
vortices, always evolves as ~ iSg, which is not consistent 
with our results. In a very recent experimental work |^ 
Anderson et al. have included a figure (Fig. 3) in which 
a two-component vortex is viewed along its axis direction 
and along an orthogonal one. Views of a vorticity free 
condensate are also included. The in-trap pictures of Ref. 
1^ cannot rule out a healing length core because such a 
feature would be too small ( ^ 0.7/im ) to resolve. So 
in this experiment Chui's condensates should look like a 
vortex free condensate. As it may be seen in the above 
mentioned figure the images of a system with and without 
vorticity are quite different. In particular, we want to 
note that in the vortex case there is a visible evidence 
of a cylinder type segregation. However to be conclusive 
about these statements more experimental information is 
required. 

In Refs. 1^,^ the authors employed, for describing the 
stationary states, wave functions expanded on a finite 
subset of the harmonic oscillator basis. This means that 
they did not work with exact solutions; unfortunately in 
these articles there is no precise information about the 
structure of the states. However, looking at both the 
domains of the wave function (Fig. 2) and the form of 
the density itself (Fig. 3), it seems that they should have 
taken into account a large number of elements of this 
basis to describe similar configurations. 

It is well known that the TFA turns out to be very 
accurate when systems with a large number of particles 
are considered, so the main features of our analysis could 
be regarded as a guide to check the accuracy of different 
approximations. 

Finally, we want to mention that very interesting ex- 
perimental |l^ and theoretical jl^ work dealing with 
off-centered vortices in two-species condensates, is forth- 
coming. 

In conclusion, we think that our analysis of states 
within the TFA could be helpful both from the exper- 
imental and theoretical viewpoints. In the first case, to 
prepare a configuration as close as possible to a station- 
ary state, and in the second, to make a suitable choice 
for the set of wave functions in which such kinds of states 
are to be expanded. 

We gratefully acknowledge useful conversations with 
B. P. Anderson. 
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FIG. 1. Domains of the different families of surfaces Sq 
in parameter space {A,B). In each region we indicate the re- 
lated z{r) curve. The thick line B — and A > corresponds 
to spheres. The B > region corresponds to cylinder-like sur- 
faces. And for the conditions A > 0, B < 0, and for B > 
the surfaces are toroids. The dashed line determines the sep- 
aration between these families of surfaces and empty sets. 
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FIG. 2. Density contours in the (r,z) plane for compo- 
nent |1) (1 2)) in the first (second) row of subfigures. The 
columns correspond to the following situations: (left) no 
vortex (7711,7712) = (0,0), (middle) a vortex in |1) species 
(7711,7712) — (1,0), and (right) a vortex in the less repulsive 
component (7711,7712) — (0, 1). The density contour spacing is 
10"^ ^m-^ 
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FIG. 3. Density profiles |*ip (solid line) and |*2|^ 
(dashed line) as function of r at z = for the same three dis- 
tributions of vorticity described in the previous figure. The 
density is given in units of 10"* fim~^. 
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